NEFSC Exploratory Data Analysis

Exporatory Data Analysis of NEFSC Bottom Trawl Survey Data

Exploration of spatial and temporal patterns in abundance, and bodymass of fishes from the Northeast groundfish survey.

# Do some formatting
trawldat <- weights_20 %>% 
  mutate(
    id = as.character(id),
    season = ifelse(season == "SPRING", "Spring", "Fall"),
    season = factor(season, levels = c("Spring", "Fall"))
  )

# Run Summary Functions
ann_means <- ss_annual_summary(trawldat)
seasonals <- ss_seasonal_summary(trawldat)

# bind them so you can facet
summs <- bind_rows(ann_means, seasonals) %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall", "Spring + Fall")))


#drop haddock to see if that changes the bump
haddock_ann  <- ss_annual_summary(filter(trawldat, comname != "haddock"))
haddock_seas <- ss_seasonal_summary(filter(trawldat, comname != "haddock"))
no_haddock <- bind_rows(haddock_ann, haddock_seas) %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall", "Spring + Fall")))

Temporal Trends

Patterns in annual abundance or CPUE of High Profile Species at annual and/or seasonal levels.

Total Biomass

Biomass of a given species from a station is a function of the adjusted number caught, their lengths, and the corresponding weights at-length using the species specific length-weight relationship coefficients. sex-specific coefficients and seasonal growth coefficients are used when available.

Total biomass in the following plots represent the total of all individual biomass estimates from all species with available growth coefficients, for each year and within each survey season.

L-W Derived

bio_1 <- ann_means %>% 
  ggplot(aes(est_year, lw_biomass_kg, color = season, linetype = "All Species")) +
    geom_line(group = 1) +
    geom_line(data = haddock_ann, 
              aes(est_year, lw_biomass_kg, color = season, linetype = "Haddock Removed")) +
    scale_y_continuous(labels = scales::comma_format()) +
    scale_color_manual(values = "black") +
    labs(x = "", y = "Total Biomass (kg)")

bio_2 <- seasonals %>% 
  ggplot(aes(est_year, lw_biomass_kg, color = season)) +
    geom_line() +
    scale_y_continuous(labels = scales::comma_format()) +
    labs(x = "", y = "Total Biomass (kg)")

bio_1 / bio_2

Gross Shipboard Weight

ship_1 <- ann_means %>% 
  ggplot(aes(est_year, fscs_biomass, color = season, linetype = "All Species")) +
    geom_line(group = 1) +
    geom_line(data = haddock_ann, 
              aes(est_year, fscs_biomass, color = season, linetype = "Haddock Removed")) +
    scale_color_manual(values = "black") +
    scale_y_continuous(labels = scales::comma_format()) +
    labs(x = "", y = "Total Biomass (kg)")

ship_2 <- seasonals %>% 
  ggplot(aes(est_year, fscs_biomass, color = season)) +
    geom_line() +
    scale_y_continuous(labels = scales::comma_format()) +
    labs(x = "", y = "Total Biomass (kg)")

ship_1 / ship_2

Area Stratified

strat_1 <- ann_means %>% 
  ggplot(aes(est_year, fscs_strat_biomass, color = season, linetype = "All Species")) +
    geom_line(group = 1) +
    geom_line(data = haddock_ann, 
              aes(est_year, fscs_strat_biomass, color = season, linetype = "Haddock Removed")) +
    scale_color_manual(values = "black") +
    scale_y_continuous(labels = scales::comma_format()) +
    labs(x = "", y = "Stratum Area Weighted Biomass")

strat_2 <- seasonals %>% 
  ggplot(aes(est_year, fscs_strat_biomass, color = season)) +
    geom_line() +
    scale_y_continuous(labels = scales::comma_format()) +
    labs(x = "", y = "Stratum Area Weighted Biomass")

strat_1 / strat_2

CPUE

Biomass per unit effort in the following figures represents the total biomass divided by the effort in number of stations. Adjustments for the 2008 vessel change are assumed to be corrected for in the data. There has been no correction for differences in stratum area in these figures.

cpue_1 <- ann_means %>% 
  ggplot(aes(est_year, lw_biomass_per_station, color = season, linetype = "All Species")) +
    geom_line(group = 1) +
    geom_line(data = haddock_ann, 
              aes(est_year, lw_biomass_per_station, color = season, linetype = "Haddock Removed")) +
    scale_color_manual(values = "black") +
    labs(x = "", y = "Adjusted Biomass per Station (kg)")

cpue_2<- seasonals %>% 
  ggplot(aes(est_year, lw_biomass_per_station, color = season)) +
    geom_line() +
    labs(x = "", y = "Adjusted Biomass per Station (kg)")

cpue_1 / cpue_2

Relative Contribution

For the following figures on relative contribution, axes for the two survey seasons are used to highlight relative differences in total biomass and cpue between the two survey seasons. A dashed line is overlayed to represent what an equal allocation would resemble. Years above the line indicate that more biomass/cpue occurred in the Fall and vice versa.

bio_comp <- seasonals %>% select(est_year,  season, lw_biomass_kg) %>% 
  pivot_wider(names_from = season, values_from = lw_biomass_kg) %>% 
  mutate(decade = floor_decade(est_year)) %>% 
  ggplot(aes(Spring, Fall, color = decade)) +
    geom_abline(slope = 1, intercept = 0, linetype = 2, color = "gray50", size = 1) +
    geom_point() +
    geom_mark_ellipse() +
    scale_y_continuous(labels = scales::comma_format()) +
    scale_x_continuous(labels = scales::comma_format()) +
    labs(x = "Spring", y = "Fall", subtitle = "Total Biomass (kg)")+
    facet_wrap(~decade) +
    theme(legend.position = "none")


cpue_comp <- seasonals %>% select(est_year,  season, lw_biomass_per_station) %>% 
  pivot_wider(names_from = season, values_from = lw_biomass_per_station) %>% 
  mutate(decade = floor_decade(est_year)) %>% 
  ggplot(aes(Spring, Fall, color = decade)) +
    geom_abline(slope = 1, intercept = 0, linetype = 2, color = "gray50", size = 1) +
    geom_point() +
    geom_mark_ellipse() +
    labs(x = "Spring", y = "Fall", subtitle = "Biomass (kg) / Standard Tow") +
    facet_wrap(~decade) +
    theme(legend.position = "none")

Biomass

bio_comp

CPUE

cpue_comp

Species Groups

Relative contribution to overall biomass by species group based on management and/or life history: Demersal, Elasmobranch, Groundfish, Pelagic

# Species Class Bins Data Prep
class_totals <- trawldat %>% 
  group_by(est_year, spec_class) %>% 
  summarise(`Total Class Biomass` = sum(sum_weight_kg)) %>% 
  ungroup() %>% 
  left_join(ann_means, by = "est_year") %>% 
  mutate(`Percent of Annual Biomass` = (`Total Class Biomass` / lw_biomass_kg) * 100 ,
         `Species Category` = spec_class,
         Year = as.character(est_year))

# Class Proportions
class_totals %>% 
  ggplot(aes(x = fct_rev(Year), y = `Percent of Annual Biomass`, fill = `Species Category`)) + 
    geom_col(color = "white", size = 0.2) +
    labs(x = "", y = "Proportion of Total Annual Biomass") +
    coord_flip() +
    scale_fill_gmri() +
    scale_y_continuous(position = "right", ) +
    guides("fill" = guide_legend(title.position = "top", title.hjust = 0.5)) +
    theme(legend.position = "bottom")

Mean Individual Biomass

Check these across years and also by season

summs %>% 
  ggplot(aes(est_year, mean_ind_bodymass * 1000, color = season)) +
    geom_line(show.legend = FALSE) +
    scale_color_manual(values = c(
      "Spring + Fall" = "black",
      "Spring" = as.character(gmri_cols("gmri blue")),
      "Fall"   = as.character(gmri_cols("orange")))) +
  #theme(legend.position = "bottom", legend.title = element_blank()) +
  labs(x = "", y = "Mean Individual Biomass (g)") +
  facet_wrap(~season, ncol = 1)

# # What do the histograms look like:
# trawldat %>%
#   select(id, est_year, area, season, length, ind_weight_kg, numlen_adj) %>%
#   uncount(numlen_adj) %>%
#   ggplot(aes(x = ind_weight_kg)) +
#   scale_y_continuous(labels = scales::comma_format()) +
#   geom_histogram()

Mean Individual Length

summs %>% 
  ggplot(aes(est_year, mean_ind_length, color = season)) +
    geom_line(show.legend = FALSE) +
    scale_color_manual(values = c(
      "Spring + Fall" = "black",
      "Spring" = as.character(gmri_cols("gmri blue")),
      "Fall"   = as.character(gmri_cols("orange")))) +
  #theme(legend.position = "bottom", legend.title = element_blank()) +
  labs(x = "", y = "Mean Individual Length (cm)", title = "Seems okay - but same method as weight...") +
  facet_wrap(~season, ncol = 1)

# # What do the histograms look like:
# trawldat %>%
#   select(id, est_year, area, season, length, ind_weight_kg, numlen_adj) %>%
#   uncount(numlen_adj) %>%
#   ggplot(aes(x = length)) +
#   scale_y_continuous(labels = scales::comma_format()) +
#   geom_histogram()

Top 5 Species by Weight

Since 2008 there has been a large increase in the amount of total biomass caught and the adjusted CPUE. As a way to look into the top species contributing to this I have pulled out the top 5 species for thbat year by total biomass caught.

# Pull top 5 species from each year, by biomass
top_5 <- trawldat %>% 
  filter(est_year >= 2004) %>% 
  group_by(est_year, comname) %>% 
  summarise(n_fish = sum(numlen_adj),
            lw_biomass = sum(sum_weight_kg)) %>% 
  ungroup() %>% 
  split(.$est_year) %>% 
  map_dfr(function(yr_group){
    yr_group  %>% 
      slice_max(order_by = lw_biomass, n = 5)
  })
  

# Plot each year, re-ordered by biomass
t5_plots <- top_5 %>% 
  split(.$est_year) %>% 
  map(function(x){
    x %>% mutate(comname = fct_drop(comname),
           comname = fct_reorder(comname, .fun = max, lw_biomass)) %>% 
    ggplot(aes(x = lw_biomass, y = comname)) +
      geom_point() +
      geom_segment(aes(yend = comname, x = 0, xend = lw_biomass)) +
      facet_wrap(~est_year) +
      scale_x_continuous(labels = scales::comma_format()) +
      labs(x = "Total Biomass (kg)", y = "")
})

t1 <- ((t5_plots$`2004` | t5_plots$`2005`) / (t5_plots$`2006` | t5_plots$`2007`))
t2 <- ((t5_plots$`2008` | t5_plots$`2009`) / (t5_plots$`2010` | t5_plots$`2011`))
t3 <- ((t5_plots$`2012` | t5_plots$`2013`) / (t5_plots$`2014` | t5_plots$`2015`))
t4 <- ((t5_plots$`2016` | t5_plots$`2017`) / (t5_plots$`2018` | t5_plots$`2019`))

2004 - 2007

t1

2008 - 2011

t2 

2012 - 2015

t3

2016 - 2019

t4

Spatial Patterns

For large regions like Georges Bank and the Gulf of Maine, what kind of patterns are we seeing.

# Load the strata
survey_strata <- read_sf(str_c(res_path, "Shapefiles/BottomTrawlStrata/BTS_Strata.shp"))  %>% 
  clean_names() %>% 
  filter(strata >= 01010 ,
         strata <= 01760,
         strata != 1310,
         strata != 1320,
         strata != 1330,
         strata != 1350,
         strata != 1410,
         strata != 1420,
         strata != 1490) 


# Key to which strata = which regions
strata_key <- list(
  "Georges Bank"          = as.character(13:23),
  "Gulf of Maine"         = as.character(24:40),
  "Southern New England"  = str_pad(as.character(1:12), width = 2, pad = "0", side = "left"),
  "Mid-Atlantic Bight"    = as.character(61:76))


# Assign Areas
survey_strata <- survey_strata %>% 
  mutate(
    strata = str_pad(strata, width = 5, pad = "0", side = "left"),
    strata_num = str_sub(strata, 3, 4),
    area = case_when(
      strata_num %in% strata_key$`Georges Bank` ~ "Georges Bank",
      strata_num %in% strata_key$`Gulf of Maine` ~ "Gulf of Maine",
      strata_num %in% strata_key$`Southern New England` ~ "Southern New England",
      strata_num %in% strata_key$`Mid-Atlantic Bight` ~ "Mid-Atlantic Bight",
    TRUE ~ "Outside Major Study Areas"
  )) %>% 
  select(finstr_id, strata, strata_num, area, a2, str2, set, stratuma, str3, geometry)

# Load new england
new_england <- ne_states("united states of america") %>% st_as_sf(crs = 4326) 
canada <- ne_states("canada") %>% st_as_sf(crs = 4326) 


# Make trawl data an sf dataset
trawl_sf <- trawldat %>% st_as_sf(coords = c("decdeg_beglon", "decdeg_beglat"), crs = 4326)

Trawl Regions

# Plot to check
ggplot() +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  geom_sf(data = survey_strata, aes(fill = area)) +
  coord_sf(xlim = c(-77, -65.5), ylim = c(34, 45.75), expand = FALSE) +
  guides(fill = guide_legend(nrow = 2)) +
  theme_bw() +
  theme(legend.position = "bottom", legend.title = element_blank())

Regional Summaries

# Just Area, all seasona
area_summs <- trawldat %>% 
  group_by(area) %>% 
  summarise(
    season = "Spring + Fall",
    lw_biomass_kg = sum(sum_weight_kg, na.rm = T),
    n_stations = n_distinct(id),
    lw_biomass_per_station = lw_biomass_kg / n_stations,
    mean_ind_bodymass = weighted.mean(ind_weight_kg, weights = numlen_adj),
    mean_ind_length = weighted.mean(length, weights = numlen_adj)
  )

# Area x Season
seas_area <- trawldat %>% 
  group_by(area, season) %>% 
  summarise(
    lw_biomass_kg = sum(sum_weight_kg, na.rm = T),
    n_stations = n_distinct(id),
    lw_biomass_per_station = lw_biomass_kg / n_stations,
    mean_ind_bodymass = weighted.mean(ind_weight_kg, weights = numlen_adj),
    mean_ind_length = weighted.mean(length, weights = numlen_adj)
  )


# Combine those two
summs_combined <- bind_rows(area_summs, seas_area) %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall", "Spring + Fall")))

summs_combined %>% 
  mutate_if(is.numeric,round, 2) %>% 
  arrange(area,season) %>% 
  knitr::kable()

area	season	lw_biomass_kg	n_stations	lw_biomass_per_station	mean_ind_bodymass	mean_ind_length
GB	Spring	4670398	3005	1554.21	0.88	40.14
GB	Fall	10430280	3129	3333.42	0.74	37.59
GB	Spring + Fall	15100678	6134	2461.80	0.81	38.79
GoM	Spring	3046226	3636	837.80	0.65	34.95
GoM	Fall	8135976	3747	2171.33	0.73	36.40
GoM	Spring + Fall	11182202	7383	1514.59	0.69	35.72
MAB	Spring	6591618	2550	2584.95	0.84	43.86
MAB	Fall	1299542	2497	520.44	0.64	25.53
MAB	Spring + Fall	7891159	5047	1563.53	0.78	38.08
SNE	Spring	4125352	2868	1438.41	0.57	36.42
SNE	Fall	2130115	2751	774.31	0.48	32.70
SNE	Spring + Fall	6255467	5619	1113.27	0.54	34.96

# Year x Area
area_summs_y <- trawldat %>% 
  group_by(est_year, area) %>% 
  summarise(
    season = "Spring + Fall",
    lw_biomass_kg = sum(sum_weight_kg, na.rm = T),
    n_stations = n_distinct(id),
    lw_biomass_per_station = lw_biomass_kg / n_stations,
    mean_ind_bodymass = weighted.mean(ind_weight_kg, weights = numlen_adj),
    mean_ind_length = weighted.mean(length, weights = numlen_adj)
  )

Total Biomass

area_summs_y %>% 
  ggplot(aes(est_year, lw_biomass_kg)) +
    geom_line() +
    facet_wrap(~area, ncol = 2) +
    scale_y_continuous(labels = scales::comma_format()) +
    labs(x = "", y = "Total Biomass (kg)")

CPUE

area_summs_y %>% 
  ggplot(aes(est_year, lw_biomass_per_station)) +
    geom_line() +
    facet_wrap(~area, ncol = 2) +
    labs(x = "", y = "Adjusted Biomass per Station (kg)")

Comparisons to Older Data

Concerns have been raised that the datasets obtained through the NEFSC are inconsistent in some areas over time. The following plots seek to identify differences in a dataset obtained in 2016 from what we currently are exploring with the 2020 dataset.

Each file was processed for size-spectra analysis using the same processing steps. This includes the same species codes, the same abundance and stratification adjustments, and the same L-W derived biomasses.

weights_16 <- read_csv(here::here("data/ss_prepped_data/survdat_2016_ss.csv"),
                       col_types = cols(),
                       guess_max = 1e5)
weights_19 <- read_csv(here::here("data/ss_prepped_data/survdat_2019_ss.csv"),
                       col_types = cols(),
                       guess_max = 1e5)
weights_20 <- read_csv(here::here("data/ss_prepped_data/survdat_2020_ss.csv"),
                       col_types = cols(),
                       guess_max = 1e5)

# run summaries
summ_16 <- ss_regional_differences(weights_16) %>% mutate(source = "2016")
summ_19 <- ss_regional_differences(weights_19) %>% mutate(source = "2019")
summ_20 <- ss_regional_differences(weights_20) %>% mutate(source = "2020")
reg_summs <- bind_rows(list(summ_16, summ_19, summ_20))

# Total Biomass
p1 <- reg_summs %>% 
  ggplot(aes(est_year, lw_biomass_kg, color = source)) +
  geom_line(show.legend = F) +
  scale_y_continuous(labels = scales::comma_format()) +
  facet_wrap(~area, ncol = 1, scales = "free") +
  labs(x = "", y = "Total Biomass \n (L-W Regressions)")

# Total Biomass - FSCS
p2 <- reg_summs %>% 
  ggplot(aes(est_year, fscs_biomass, color = source)) +
  geom_line() +
  scale_y_continuous(labels = scales::comma_format()) +
  facet_wrap(~area, ncol = 1) +
  labs(x = "", y = "Total Biomass \n (FSCS Haul Weights)")

# effort
p3 <- reg_summs %>% 
  ggplot(aes(est_year, n_stations, color = source)) +
  geom_line(show.legend = F) +
  scale_y_continuous(labels = scales::comma_format()) +
  facet_wrap(~area, ncol = 1) +
  labs(x = "", y = "Effort (haul count)")

# Species 
p4 <- reg_summs %>% 
  ggplot(aes(est_year, n_species, color = source)) +
  geom_line(show.legend = F) +
  scale_y_continuous(labels = scales::comma_format()) +
  facet_wrap(~area, ncol = 1) +
  labs(x = "", y = "Distinct Species")

p1 + p2 + p3 + p4

 

---

A work by Adam A. Kemberling

Akemberling@gmri.org